home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / calctree.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  21.4 KB  |  1,028 lines  |  [TEXT/R*ch]
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5. #include <stdarg.h>
  6. #include <ctype.h>
  7. #include "clustalw.h"
  8.  
  9. /*
  10.  *   Prototypes
  11.  */
  12. extern void *ckalloc(size_t bytes);
  13. extern void ckfree(void *);
  14.  
  15. void calc_seq_weights(int first_seq, int last_seq);
  16. void create_sets(int first_seq, int last_seq);
  17. int calc_similarities(int nseqs);
  18. int calc_distances(int nseqs);
  19. int read_tree(char *treefile, int first_seq, int last_seq);
  20.  
  21. void skip_space(FILE *fd);
  22. treeptr avail(void);
  23. void clear_tree(treeptr p);
  24. void set_info(treeptr p, treeptr parent, int pleaf, char *pname, float pdist);
  25. treeptr reroot(treeptr ptree);
  26. treeptr insert_root(treeptr p, float diff);
  27. float calc_root_mean(treeptr root, float *maxdist);
  28. float calc_mean(treeptr root, treeptr nptr, float *maxdist);
  29. void create_tree(treeptr ptree, treeptr parent);
  30. void create_node(treeptr pptr, treeptr parent);
  31. treeptr insert_node(treeptr pptr);
  32. void order_nodes(void);
  33. int calc_weight(int leaf);
  34. void group_seqs(treeptr root, int *groups, int nseqs);
  35. void save_set(int n, int *groups);
  36. void mark_group1(treeptr p, int *groups, int n);
  37. void mark_group2(treeptr p, int *groups, int n);
  38.  
  39. /*
  40.  *   Global variables
  41.  */
  42.  
  43. extern Boolean distance_tree;
  44. extern int debug;
  45. extern double **tmat;
  46. extern int **sets;
  47. extern int nsets;
  48. extern char **names;
  49. extern int *seq_weight;
  50.  
  51. char ch;
  52. FILE *fd;
  53. treeptr lptr[MAXN];
  54. treeptr olptr[MAXN];
  55. treeptr nptr[MAXTREE];
  56. treeptr ptrs[MAXTREE];
  57. int nnodes = 0;
  58. int ntotal = 0;
  59. int rooted_tree = TRUE;
  60. static treeptr seq_tree,root;
  61. static int *groups, numseq;
  62.  
  63. void calc_seq_weights(int first_seq, int last_seq)
  64. {
  65.   int   i, j, k, nseqs;
  66.   int   temp, sum, *weight;
  67.  
  68.  
  69. /*
  70.   If there are more than three sequences....
  71. */
  72.   nseqs = last_seq-first_seq;
  73.   if ((nseqs > 3) && (distance_tree == TRUE))
  74.      {
  75. /*
  76.   Calculate sequence weights based on Phylip tree.
  77. */
  78.       weight = (int *)ckalloc((nseqs+1) * sizeof(int));
  79.  
  80.       for (i=first_seq; i<last_seq; i++)
  81.            weight[i] = calc_weight(i);
  82.  
  83. /*
  84.   Normalise the weights, such that the sum of the weights = 1000
  85. */
  86.  
  87.          sum = 0;
  88.          for (i=first_seq; i<last_seq; i++)
  89.             sum += weight[i];
  90.  
  91.          if (sum == 0)
  92.           {
  93.             for (i=first_seq; i<last_seq; i++)
  94.                weight[i] = 1;
  95.             sum = i;
  96.           }
  97.  
  98.          for (i=first_seq; i<last_seq; i++)
  99.            {
  100.               seq_weight[i] = (weight[i] * 1000) / sum;
  101.               if (seq_weight[i] < 1) seq_weight[i] = 1;
  102.            }
  103.  
  104.        ckfree((void *)weight);
  105.  
  106.      }
  107.  
  108.    else
  109.      {
  110. /*
  111.   Otherwise, use identity weights.
  112. */
  113.         temp = 1000 / nseqs;
  114.         for (i=first_seq; i<last_seq; i++)
  115.            seq_weight[i] = temp;
  116.      }
  117.  
  118. }
  119.  
  120. void create_sets(int first_seq, int last_seq)
  121. {
  122.   int   i, j, k, nseqs;
  123.  
  124.   nsets = 0;
  125.   nseqs = last_seq-first_seq;
  126.   if (nseqs > 3)
  127.      {
  128. /*
  129.   If there are more than three sequences....
  130. */
  131.        groups = (int *)ckalloc((nseqs+1) * sizeof(int));
  132.        group_seqs(root, groups, nseqs);
  133.        ckfree((void *)groups);
  134.  
  135.      }
  136.  
  137.    else
  138.      {
  139.        groups = (int *)ckalloc((nseqs+1) * sizeof(int));
  140.        for (i=0;i<nseqs-1;i++)
  141.          {
  142.            for (j=0;j<nseqs;j++)
  143.               if (j<=i) groups[j] = 1;
  144.               else if (j==i+1) groups[j] = 2;
  145.               else groups[j] = 0;
  146.            save_set(nseqs, groups);
  147.          }
  148.        ckfree((void *)groups);
  149.      }
  150.  
  151. }
  152.  
  153. int read_tree(char *treefile, int first_seq, int last_seq)
  154. {
  155.  
  156.   char c;
  157.   char name1[MAXNAMES+1], name2[MAXNAMES+1];
  158.   int i, j, k, n, value, rootnode, rootdirection;
  159.   int found;
  160.  
  161.   numseq = 0;
  162.   nnodes = 0;
  163.   ntotal = 0;
  164.   rooted_tree = TRUE;
  165.  
  166. #ifdef VMS
  167.   if ((fd = fopen(treefile,"r","rat=cr","rfm=var")) == NULL)
  168. #else
  169.   if ((fd = fopen(treefile, "r")) == NULL)
  170. #endif
  171.     {
  172.       fprintf(stderr, "Error: cannot open %s\n", treefile);
  173.       return(0);
  174.     }
  175.  
  176.   skip_space(fd);
  177.   ch = (char)getc(fd);
  178.   if (ch != '(')
  179.     {
  180.       fprintf(stderr, "Error: Wrong format in tree file %s\n", treefile);
  181.       return(0);
  182.     }
  183.   rewind(fd);
  184.  
  185.   distance_tree = TRUE;
  186.  
  187.   seq_tree = avail();
  188.   set_info(seq_tree, NULL, 0, "", 0.0);
  189.  
  190.   create_tree(seq_tree,NULL);
  191.   fclose(fd);
  192.  
  193. /*
  194.   If the tree is unrooted, reroot the tree - ie. minimise the difference
  195.   between the mean root->leaf distances for the left and right branches of
  196.   the tree.
  197. */
  198.  
  199.   if (distance_tree == FALSE)
  200.      {
  201.       if (rooted_tree == FALSE)
  202.           {
  203.                 fprintf(stderr,
  204.              "Error: input tree is unrooted and has no distances, cannot align sequences\n");
  205.              return(0);
  206.           }
  207.      }
  208.  
  209.   if (rooted_tree == FALSE)
  210.      {
  211.         root = reroot(seq_tree);
  212.      }
  213.   else
  214.      {
  215.         root = seq_tree;
  216.      }
  217.  
  218. /*
  219.   calculate the 'order' of each node.
  220. */
  221.   order_nodes();
  222.  
  223.  
  224.   if (numseq != last_seq-first_seq)
  225.      {
  226.          fprintf(stderr," Error: tree not compatible with alignment (%d sequences in alignment and %d in tree\n",
  227.          (pint)last_seq-first_seq,(pint)numseq);
  228.          return(0);
  229.      }
  230.  
  231.   if (numseq > 3)
  232.      {
  233. /*
  234.   If there are more than three sequences....
  235. */
  236. /*
  237.   assign the sequence nodes (in the same order as in the alignment file)
  238. */
  239.       for (i=first_seq; i<last_seq; i++)
  240.        {
  241.          if (strlen(names[i+1]) > MAXNAMES)
  242.              fprintf(stderr,"Warning: name %s is too long for PHYLIP tree format (max 10 chars)\n", names[i+1]);
  243.  
  244.          for (k=0; k< strlen(names[i+1]) && k<MAXNAMES ; k++)
  245.            {
  246.              c = names[i+1][k];
  247.              if ((c>0x40) && (c<0x5b)) c=c | 0x20;
  248.              if (c == ' ') c = '_';
  249.              name2[k] = c;
  250.            }
  251.          name2[k]='\0';
  252.          found = FALSE;
  253.          for (j=0; j<numseq; j++)
  254.            {
  255.             for (k=0; k< strlen(lptr[j]->name) && k<MAXNAMES ; k++)
  256.               {
  257.                 c = lptr[j]->name[k];
  258.                 if ((c>0x40) && (c<0x5b)) c=c | 0x20;
  259.                 name1[k] = c;
  260.               }
  261.             name1[k]='\0';
  262.             if (strcmp(name1, name2) == 0)
  263.               {
  264.                 olptr[i] = lptr[j];
  265.                 found = TRUE;
  266.               }
  267.            }
  268.          if (found == FALSE)
  269.            {
  270.              fprintf(stderr," Error: tree not compatible with alignment: %s not found\n",
  271.                                       name2);
  272.              return(0);
  273.            }
  274.        }
  275.  
  276.      }
  277.    return(1);
  278. }
  279.  
  280. void create_tree(treeptr ptree, treeptr parent)
  281. {
  282.    treeptr p;
  283.  
  284.    int i, type;
  285.    float dist;
  286.    char name[MAXNAMES+1];
  287.  
  288. /*
  289.   is this a node or a leaf ?
  290. */
  291.   skip_space(fd);
  292.   ch = (char)getc(fd);
  293.   if (ch == '(')
  294.     {  
  295. /*
  296.    this must be a node....
  297. */
  298.       type = NODE;
  299.       name[0] = '\0';
  300.       ptrs[ntotal] = nptr[nnodes] = ptree;
  301.       nnodes++;
  302.       ntotal++;
  303.  
  304.       create_node(ptree, parent);
  305.  
  306.       p = ptree->left;
  307.       create_tree(p, ptree);
  308.            
  309.       if ( ch == ',')
  310.        {
  311.           p = ptree->right;
  312.           create_tree(p, ptree);
  313.           if ( ch == ',')
  314.             {
  315.                ptree = insert_node(ptree);
  316.                ptrs[ntotal] = nptr[nnodes] = ptree;
  317.                nnodes++;
  318.                ntotal++;
  319.                p = ptree->right;
  320.                create_tree(p, ptree);
  321.                rooted_tree = FALSE;
  322.             }
  323.        }
  324.  
  325.       skip_space(fd);
  326.       ch = (char)getc(fd);
  327.     }
  328. /*
  329.    ...otherwise, this is a leaf
  330. */
  331.   else
  332.     {
  333.       type = LEAF;
  334.       ptrs[ntotal++] = lptr[numseq++] = ptree;
  335. /*
  336.    get the sequence name
  337. */
  338.       name[0] = ch;
  339.       ch = (char)getc(fd);
  340.       i = 1;
  341.       while ((ch != ':') && (ch != ',') && (ch != ')'))
  342.         {
  343.           if (i < MAXNAMES) name[i++] = ch;
  344.           ch = (char)getc(fd);
  345.         }
  346.       name[i] = '\0';
  347.       if (ch != ':')
  348.          {
  349.            distance_tree = FALSE;
  350.            dist = 0.0;
  351.          }
  352.     }
  353.  
  354. /*
  355.    get the distance information
  356. */
  357.   dist = 0.0;
  358.   if (ch == ':')
  359.      {
  360.        skip_space(fd);
  361.        fscanf(fd,"%f",&dist);
  362.        skip_space(fd);
  363.        ch = (char)getc(fd);
  364.      }
  365.    set_info(ptree, parent, type, name, dist);
  366.  
  367.  
  368. }
  369.  
  370. void create_node(treeptr pptr, treeptr parent)
  371. {
  372.   treeptr t;
  373.  
  374.   pptr->parent = parent;
  375.   t = avail();
  376.   pptr->left = t;
  377.   t = avail();
  378.   pptr->right = t;
  379.     
  380. }
  381.  
  382. treeptr insert_node(treeptr pptr)
  383. {
  384.  
  385.    treeptr newnode;
  386.  
  387.    newnode = avail();
  388.    create_node(newnode, pptr->parent);
  389.  
  390.    newnode->left = pptr;
  391.    pptr->parent = newnode;
  392.  
  393.    set_info(newnode, pptr->parent, NODE, "", 0.0);
  394.  
  395.    return(newnode);
  396. }
  397.  
  398. void skip_space(FILE *fd)
  399. {
  400.   int   c;
  401.   
  402.   do
  403.      c = getc(fd);
  404.   while(isspace(c));
  405.  
  406.   ungetc(c, fd);
  407. }
  408.  
  409. treeptr avail(void)
  410. {
  411.    treeptr p;
  412.    p = ckalloc(sizeof(stree));
  413.    p->left = NULL;
  414.    p->right = NULL;
  415.    p->parent = NULL;
  416.    p->dist = 0.0;
  417.    p->leaf = 0;
  418.    p->order = 0;
  419.    p->name[0] = '\0';
  420.    return(p);
  421. }
  422.  
  423. void clear_tree(treeptr p)
  424. {
  425.    if (p==NULL) p = root;
  426.    if (p->left != NULL)
  427.      {
  428.        clear_tree(p->left);
  429.      }
  430.    if (p->right != NULL)
  431.      {
  432.        clear_tree(p->right);
  433.      }
  434.    p->left = NULL;
  435.    p->right = NULL;
  436.    ckfree((void *)p);
  437. }
  438.  
  439. void set_info(treeptr p, treeptr parent, int pleaf, char *pname, float pdist)
  440. {
  441.    p->parent = parent;
  442.    p->leaf = pleaf;
  443.    p->dist = pdist;
  444.    p->order = 0;
  445.    strcpy(p->name, pname);
  446.    if (p->leaf == TRUE)
  447.      {
  448.         p->left = NULL;
  449.         p->right = NULL;
  450.      }
  451. }
  452.  
  453. treeptr reroot(treeptr ptree)
  454. {
  455.  
  456.    treeptr p,q,t, rootnode, rootptr;
  457.    float   lmean, rmean;
  458.    float   dist, diff, mindiff, mindepth = 1.0, maxdist;
  459.    int     i,j;
  460.    int     first = TRUE;
  461.  
  462. /*
  463.   find the difference between the means of leaf->node
  464.   distances on the left and on the right of each node
  465. */
  466.    rootptr = ptree;
  467.    for (i=0; i<ntotal; i++)
  468.      {
  469.         p = ptrs[i];
  470.         if (p->parent == NULL)
  471.            diff = calc_root_mean(p, &maxdist);
  472.         else
  473.            diff = calc_mean(ptree, p, &maxdist);
  474.  
  475.         if ((diff == 0) || ((diff > 0) && (diff < 2 * p->dist)))
  476.           {
  477.               if ((maxdist < mindepth) || (first == TRUE))
  478.                  {
  479.                     first = FALSE;
  480.                     rootptr = p;
  481.                     mindepth = maxdist;
  482.                     mindiff = diff;
  483.                  }
  484.            }
  485.  
  486.      }
  487.  
  488. /*
  489.   insert a new node as the ancestor of the node which produces the shallowest
  490.   tree.
  491. */
  492.    if (rootptr == ptree)
  493.      {
  494.         mindiff = rootptr->left->dist + rootptr->right->dist;
  495.         rootptr = rootptr->right;
  496.      }
  497.    rootnode = insert_root(rootptr, mindiff);
  498.   
  499.    diff = calc_root_mean(rootnode, &maxdist);
  500.  
  501.    return(rootnode);
  502. }
  503.  
  504. treeptr insert_root(treeptr p, float diff)
  505. {
  506.    treeptr newp, prev, prevp, q, t;
  507.    float dist, prevdist, prevdiff, td;
  508.  
  509.    newp = avail();
  510.  
  511.    t = p->parent;
  512.    prevdist = t->dist;
  513.  
  514.    p->parent = newp;
  515.  
  516.    dist = p->dist;
  517.  
  518.    p->dist = diff / 2;
  519.    if (p->dist < 0.0) p->dist = 0.0;
  520.    if (p->dist > dist) p->dist = dist;
  521.  
  522.    t->dist = dist - p->dist; 
  523.  
  524.    newp->left = t;
  525.    newp->right = p;
  526.    newp->parent = NULL;
  527.    newp->dist = 0.0;
  528.    newp->leaf = NODE;
  529.  
  530.    if (t->left == p) t->left = t->parent;
  531.    else t->right = t->parent;
  532.  
  533.    prev = t;
  534.    q = t->parent;
  535.  
  536.    t->parent = newp;
  537.  
  538.    while (q != NULL)
  539.      {
  540.         if (q->left == prev)
  541.            {
  542.               q->left = q->parent;
  543.               q->parent = prev;
  544.               td = q->dist;
  545.               q->dist = prevdist;
  546.               prevdist = td;
  547.               prev = q;
  548.               q = q->left;
  549.            }
  550.         else
  551.            {
  552.               q->right = q->parent;
  553.               q->parent = prev;
  554.               td = q->dist;
  555.               q->dist = prevdist;
  556.               prevdist = td;
  557.               prev = q;
  558.               q = q->right;
  559.            }
  560.     }
  561.  
  562. /*
  563.    remove the old root node
  564. */
  565.    q = prev;
  566.    if (q->left == NULL)
  567.       {
  568.          dist = q->dist;
  569.          q = q->right;
  570.          q->dist += dist;
  571.          q->parent = prev->parent;
  572.          if (prev->parent->left == prev)
  573.             prev->parent->left = q;
  574.          else
  575.             prev->parent->right = q;
  576.          prev->right = NULL;
  577.       }
  578.    else
  579.       {
  580.          dist = q->dist;
  581.          q = q->left;
  582.          q->dist += dist;
  583.          q->parent = prev->parent;
  584.          if (prev->parent->left == prev)
  585.             prev->parent->left = q;
  586.          else
  587.             prev->parent->right = q;
  588.          prev->left = NULL;
  589.       }
  590.  
  591.    return(newp);
  592. }
  593.  
  594. float calc_root_mean(treeptr root, float *maxdist)
  595. {
  596.    float dist , lsum = 0.0, rsum = 0.0, lmean,rmean,diff;
  597.    treeptr p;
  598.    int depth = 0, i,j , n;
  599.    int nl, nr;
  600.    int direction, found;
  601. /*
  602.    for each leaf, determine whether the leaf is left or right of the root.
  603. */
  604.    dist = (*maxdist) = 0;
  605.    nl = nr = 0;
  606.    for (i=0; i< numseq; i++)
  607.      {
  608.          p = lptr[i];
  609.          dist = 0.0;
  610.          while (p->parent != root)
  611.            {
  612.                dist += p->dist;
  613.                p = p->parent;
  614.            }
  615.          if (p == root->left) direction = LEFT;
  616.          else direction = RIGHT;
  617.          dist += p->dist;
  618.  
  619.          if (direction == LEFT)
  620.            {
  621.              lsum += dist;
  622.              nl++;
  623.            }
  624.          else
  625.            {
  626.              rsum += dist;
  627.              nr++;
  628.            }
  629.         if (dist > (*maxdist)) *maxdist = dist;
  630.      }
  631.  
  632.    lmean = lsum / nl;
  633.    rmean = rsum / nr;
  634.  
  635.    diff = lmean - rmean;
  636.    return(diff);
  637. }
  638.  
  639.  
  640. float calc_mean(treeptr root, treeptr nptr, float *maxdist)
  641. {
  642.    float dist , lsum = 0.0, rsum = 0.0, lmean,rmean,diff;
  643.    treeptr p, path2root[MAXN];
  644.    float dist2node[MAXN];
  645.    int depth = 0, i,j , n;
  646.    int nl , nr;
  647.    int direction, found;
  648.  
  649. /*
  650.    determine all nodes between the selected node and the root;
  651. */
  652.    depth = (*maxdist) = dist = 0;
  653.    nl = nr = 0;
  654.    p = nptr;
  655.    while (p != NULL)
  656.      {
  657.          path2root[depth] = p;
  658.          dist += p->dist;
  659.          dist2node[depth] = dist;
  660.          p = p->parent;
  661.          depth++;
  662.      }
  663.  
  664. /*
  665.    *nl = *nr = 0;
  666.    for each leaf, determine whether the leaf is left or right of the node.
  667.    (RIGHT = descendant, LEFT = not descendant)
  668. */
  669.    for (i=0; i< numseq; i++)
  670.      {
  671.        p = lptr[i];
  672.        if (p == nptr)
  673.          {
  674.             direction = RIGHT;
  675.             dist = 0.0;
  676.          }
  677.        else
  678.          {
  679.          direction = LEFT;
  680.          dist = 0.0;
  681. /*
  682.    find the common ancestor.
  683. */
  684.          found = FALSE;
  685.          n = 0;
  686.          while ((found == FALSE) && (p->parent != NULL))
  687.            {
  688.                for (j=0; j< depth; j++)
  689.                  if (p->parent == path2root[j])
  690.                     { 
  691.                       found = TRUE;
  692.                       n = j;
  693.                     }
  694.                dist += p->dist;
  695.                p = p->parent;
  696.            }
  697.          if (p == nptr) direction = RIGHT;
  698.          }
  699.  
  700.          if (direction == LEFT)
  701.            {
  702.              lsum += dist;
  703.              lsum += dist2node[n-1];
  704.              nl++;
  705.            }
  706.          else
  707.            {
  708.              rsum += dist;
  709.              nr++;
  710.            }
  711.  
  712.         if (dist > (*maxdist)) *maxdist = dist;
  713.      }
  714.  
  715.    lmean = lsum / nl;
  716.    rmean = rsum / nr;
  717.    
  718.    diff = lmean - rmean;
  719.    return(diff);
  720. }
  721.  
  722. void order_nodes(void)
  723. {
  724.    int i;
  725.    treeptr p;
  726.  
  727.    for (i=0; i<numseq; i++)
  728.      {
  729.         p = lptr[i];
  730.         while (p != NULL)
  731.           {
  732.              p->order++;
  733.              p = p->parent;
  734.           }
  735.      }
  736. }
  737.  
  738.  
  739. int calc_weight(int leaf)
  740. {
  741.  
  742.   treeptr p;
  743.   float weight = 0.0;
  744.  
  745.   p = olptr[leaf];
  746.   while (p->parent != NULL)
  747.     {
  748.        weight += p->dist / p->order;
  749.        p = p->parent;
  750.     }
  751.  
  752.   weight *= 100.0;
  753.  
  754.   return((int)weight);
  755.  
  756. }
  757.  
  758. void group_seqs(treeptr p, int *next_groups, int nseqs)
  759. {
  760.     int i, n;
  761.     int *tmp_groups;
  762.  
  763.     tmp_groups = (int *)ckalloc((nseqs+1) * sizeof(int));
  764.     for (i=0;i<nseqs;i++)
  765.          tmp_groups[i] = 0;
  766.  
  767.     if (p->left != NULL)
  768.       {
  769.          if (p->left->leaf == NODE)
  770.             {
  771.                group_seqs(p->left, next_groups, nseqs);
  772.                for (i=0;i<nseqs;i++)
  773.                  if (next_groups[i] != 0) tmp_groups[i] = 1;
  774.             }
  775.          else
  776.             {
  777.                mark_group1(p->left, tmp_groups, nseqs);
  778.             }
  779.                
  780.       }
  781.  
  782.     if (p->right != NULL)
  783.       {
  784.          if (p->right->leaf == NODE)
  785.             {
  786.                group_seqs(p->right, next_groups, nseqs);
  787.                for (i=0;i<nseqs;i++)
  788.                     if (next_groups[i] != 0) tmp_groups[i] = 2;
  789.             }
  790.          else 
  791.             {
  792.                mark_group2(p->right, tmp_groups, nseqs);
  793.             }
  794.          save_set(nseqs, tmp_groups);
  795.       }
  796.     for (i=0;i<nseqs;i++)
  797.       next_groups[i] = tmp_groups[i];
  798.  
  799.     ckfree((void *)tmp_groups);
  800.  
  801. }
  802.  
  803. void mark_group1(treeptr p, int *groups, int n)
  804. {
  805.     int i;
  806.  
  807.     for (i=0;i<n;i++)
  808.        {
  809.          if (olptr[i] == p)
  810.               groups[i] = 1;
  811.          else
  812.               groups[i] = 0;
  813.        }
  814. }
  815.  
  816. void mark_group2(treeptr p, int *groups, int n)
  817. {
  818.     int i;
  819.  
  820.     for (i=0;i<n;i++)
  821.        {
  822.          if (olptr[i] == p)
  823.               groups[i] = 2;
  824.          else if (groups[i] != 0)
  825.               groups[i] = 1;
  826.        }
  827. }
  828.  
  829. void save_set(int n, int *groups)
  830. {
  831.     int i;
  832.  
  833.     for (i=0;i<n;i++)
  834.       sets[nsets+1][i+1] = groups[i];
  835.     nsets++;
  836. }
  837.  
  838.  
  839.  
  840. int calc_similarities(int nseqs)
  841. {
  842.    int depth = 0, i,j, k, n;
  843.    int direction, found;
  844.    treeptr p, *path2root;
  845.    float dist , sum = 0.0;
  846.    float *dist2node;
  847.    double **dmat;
  848.  
  849.    path2root = (treeptr *)ckalloc((nseqs) * sizeof(treeptr));
  850.    dist2node = (float *)ckalloc((nseqs) * sizeof(float));
  851.    dmat = (double **)ckalloc((nseqs) * sizeof(double *));
  852.    for (i=0;i<nseqs;i++)
  853.      dmat[i] = (double *)ckalloc((nseqs) * sizeof(double));
  854.  
  855.    if (nseqs > 3)
  856.     {
  857. /*
  858.    for each leaf, determine all nodes between the leaf and the root;
  859. */
  860.       for (i = 0;i<nseqs; i++)
  861.        { 
  862.           depth = dist = 0;
  863.           p = olptr[i];
  864.           while (p != NULL)
  865.             {
  866.                 path2root[depth] = p;
  867.                 dist += p->dist;
  868.                 dist2node[depth] = dist;
  869.                 p = p->parent;
  870.                 depth++;
  871.             }
  872.  
  873. /*
  874.    for each pair....
  875. */
  876.           for (j=0; j < i; j++)
  877.             {
  878.               p = olptr[j];
  879.               dist = 0.0;
  880. /*
  881.    find the common ancestor.
  882. */
  883.               found = FALSE;
  884.               n = 0;
  885.               while ((found == FALSE) && (p->parent != NULL))
  886.                 {
  887.                     for (k=0; k< depth; k++)
  888.                       if (p->parent == path2root[k])
  889.                          { 
  890.                            found = TRUE;
  891.                            n = k;
  892.                          }
  893.                     dist += p->dist;
  894.                     p = p->parent;
  895.                 }
  896.    
  897.               dmat[i][j] = dist + dist2node[n-1];
  898.             }
  899.         }
  900.  
  901.         dist = 0.0;
  902.         for (i=0;i<nseqs;i++)
  903.           {
  904.              dmat[i][i] = 0.0;
  905.              for (j=0;j<i;j++)
  906.                {
  907.                   if (dmat[i][j] < 0.01) dmat[i][j] = 0.01;
  908.                   if ((dmat[i][j] > 1.0) && (dmat[i][j] < 1.1))
  909.                           dmat[i][j] = 1.0;
  910.                   if (dmat[i][j] > dist) dist = dmat[i][j];
  911.                }
  912.           }
  913.         if (dist > 1.01) 
  914.           {
  915.              fprintf(stderr,"Warning: distances in tree are out of range (all distances must be between 0.0 and 1.0 (biggest dist is %f))\n",dist);
  916.              return(0);
  917.           }
  918.      }
  919.    else
  920.      {
  921.         for (i=0;i<nseqs;i++)
  922.           {
  923.              for (j=0;j<i;j++)
  924.                {
  925.                   dmat[i][j] = tmat[i+1][j+1];
  926.                }
  927.           }
  928.      }
  929.  
  930.    ckfree((void *)path2root);
  931.    ckfree((void *)dist2node);
  932.    for (i=0;i<nseqs;i++)
  933.      {
  934.         tmat[i][i] = 0.0;
  935.         for (j=0;j<i;j++)
  936.           {
  937.              tmat[i+1][j+1] = 100.0 - (dmat[i][j]) * 100.0;
  938.              tmat[j+1][i+1] = tmat[i+1][j+1];
  939.           }
  940.      }
  941.  
  942.    for (i=0;i<nseqs;i++) ckfree((void *)dmat[i]);
  943.    ckfree((void *)dmat);
  944.  
  945.    return(1);
  946. }
  947.  
  948. int calc_distances(int nseqs)
  949. {
  950.    int depth = 0, i,j, k, n;
  951.    int direction, found;
  952.    treeptr p, *path2root;
  953.    float dist , sum = 0.0;
  954.    float *dist2node;
  955.  
  956.    path2root = (treeptr *)ckalloc((nseqs) * sizeof(treeptr));
  957.    dist2node = (float *)ckalloc((nseqs) * sizeof(float));
  958.  
  959.    if (nseqs > 3)
  960.     {
  961. /*
  962.    for each leaf, determine all nodes between the leaf and the root;
  963. */
  964.       for (i = 0;i<nseqs; i++)
  965.        { 
  966.           depth = dist = 0;
  967.           p = olptr[i];
  968.           while (p != NULL)
  969.             {
  970.                 path2root[depth] = p;
  971.                 dist += p->dist;
  972.                 dist2node[depth] = dist;
  973.                 p = p->parent;
  974.                 depth++;
  975.             }
  976.  
  977. /*
  978.    for each pair....
  979. */
  980.           for (j=0; j < i; j++)
  981.             {
  982.               p = olptr[j];
  983.               dist = 0.0;
  984. /*
  985.    find the common ancestor.
  986. */
  987.               found = FALSE;
  988.               n = 0;
  989.               while ((found == FALSE) && (p->parent != NULL))
  990.                 {
  991.                     for (k=0; k< depth; k++)
  992.                       if (p->parent == path2root[k])
  993.                          { 
  994.                            found = TRUE;
  995.                            n = k;
  996.                          }
  997.                     dist += p->dist;
  998.                     p = p->parent;
  999.                 }
  1000.    
  1001.               tmat[i+1][j+1] = dist + dist2node[n-1];
  1002.             }
  1003.         }
  1004.  
  1005.         dist = 0.0;
  1006.         for (i=0;i<nseqs;i++)
  1007.           {
  1008.              tmat[i+1][i+1] = 0.0;
  1009.              for (j=0;j<i;j++)
  1010.                {
  1011.                   if (tmat[i+1][j+1] < 0.01) tmat[i+1][j+1] = 0.01;
  1012.                   if (tmat[i+1][j+1] > dist) dist = tmat[i+1][j+1];
  1013.                   tmat[j+1][i+1] = tmat[i+1][j+1];
  1014.                }
  1015.           }
  1016.         if (dist > 1.0) 
  1017.           {
  1018.              fprintf(stderr,"Warning: distances in tree are out of range (all distances must be between 0.0 and 1.0)\n");
  1019.              return(0);
  1020.           }
  1021.      }
  1022.  
  1023.    ckfree((void *)path2root);
  1024.    ckfree((void *)dist2node);
  1025.  
  1026.    return(1);
  1027. }
  1028.